number.of.neutral.patterns <- NULL
sw <- 0
# Sweep div prob
for (divprob in 0:4){
    # Sweep selective coefficient
    for (selcoeff in 0:5){
	# Sweep selective mut rate
	for (selmu in 0:3){
	    # Sweep neutral mut rate
	    for (neutmu in 0:2){
		# Sweep replicates
		for (replicates in 0:2){
		    sw <- sw + 1
		    # Load state of grid at end of simulation
		    grid.log <- read.csv(paste("../sim",sw,"/sim",sw,".grid.log",sep=""),header=FALSE)		    
		    grid.log <- grid.log[grid.log[,1]==7300,]
		    # -1 for NA's corresponding to dead crypts
		    number.of.neutral.patterns[sw] <- length(unique(grid.log[,5])-1)
		    print (paste("Sim",sw,"done."))
		}
	    }
	}	
    }
}

write.table(number.of.neutral.patterns,"simulation.number.and.number.of.neutral.patterns.csv",quote=FALSE,sep=",",col.names=FALSE,row.names=FALSE)


library(gplots)

# Read no of neut pat per sim
number.of.neutral.patterns <- read.csv("simulation.number.and.number.of.neutral.patterns.csv",header=FALSE)

divproblab <- c("r=0","r=0.25","r=0.5","r=0.75","r=1")
selcoefflab <- c("s=0.001","s=0.01","s=0.1","s=0.5","s=1","s=2")
selmulab <- c("10^-5","10^-6","10^-7","10^-8")
neutmulab <- c("10^-4","10^-5","10^-7")

res <- matrix(nrow=72,ncol=5)
sw <- 0
for (divprob in 1:5){
for (selcoeff in 1:6){
rrow <- 0
for (selmu in 1:4){
for (neutmu in 1:3){
#for (reps in 1:3){
    sw <- sw + 1
    res[rrow+selcoeff,divprob] <- mean(c(number.of.neutral.patterns[sw,1],number.of.neutral.patterns[sw+1,1],number.of.neutral.patterns[sw+2,1]))
    rrow <- rrow + 6
    sw <- sw + 2
#}
}
}
}
}

res <- round(res)


pdf("number.of.neutral.patterns.pdf",width=6, height=6)
sw <- 0
for (selmu in 1:4){
for (neutmu in 1:3){
sw <- sw + 1
subres <- res[(sw*6-5):(sw*6),] 
colnames(subres) <- c("r=0","r=0.25","r=0.5","r=0.75","r=1")
rownames(subres) <- c("s=0.001","s=0.01","s=0.1","s=0.5","s=1","s=2")

mp <- barplot2(subres, beside = TRUE,
        col = gray(1-(0:5/5)),        
        main = paste("smu =", selmulab[selmu],"nmu =",neutmulab[neutmu]), font.main = 4,
        plot.grid = TRUE)
#legend(0,max(subres)*9/10,legend=rownames(subres),fill=gray(1-(0:5/5)))
mtext(side=2, text="Average number of patterns per 3 runs",line=2)
mtext(side=1, text="Crypt division probability (r)",line=2)
}
}
dev.off()

